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Abstract 


The steady-state behavior of a particular type of digital phase-locked 
loop (DPLL) with integrate-and-dump circuit following the phase detector is 
characterized in terms of the probability density function (pdf) of the phase 
error in the loop. Although the loop is entirely digital from an implementa- 
tion standpoint, it operates at two extremely different sampling rates. In 
particular, the combination of a phase detector and an integrate-and-dump 
circuit operates at a very high rate whereas the loop update rate is very slow 
by comparison. Because of this dichotomy, the loop can be analyzed by hybrid 
analog/digital (s/z domain) techniques. The loop is modeled in such a general 
fashion that previous analyses of the Real-Time Combiner (RTC), Subcarrier 
Demodulator Assembly (SDA), and Symbol Synchronization Assembly (SSA) fall out 
as special cases. 
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I. 


INTRODUCTION 


The purpose of this analysis is to derive the steady state probability 
density function (pdf) of the phase error in a digital phase-locked loop 
(DPLL). The loop considered in this analysis has an integrate-and-dump device 
after the phase detector. This device together with the phase detector 
operates at a much higher rate than the rest of the loop, which makes it 
possible (or imperative) to model this part of the digital loop as an analog 
device. In previous analyses. Refs. 1, 2, and 3, each component of the loop 
was modeled either in the s or the z domain, which gave a hybrid s/z loop 
model. The resulting open and closed loop transfer functions related 
integrated phase error to the integrated input phase process. 

In the present analysis, difference equations are directly written for 
the different components of the loop. These difference equations are combined 
in order to have an overall mathematical model of the loop in the presence of 
noise. It is shown that essentially the same open/closed loop transfer func- 
tions are obtained as in the previous analyses. The main difference is that 
now the transfer functions relate the phase error to the input phase process. 
This fact permits the derivation of the steady-state pdf of the phase error. 

A method is given to transform high-order nonlinear stochastic difference 
equations. It is shown how this analysis applies to the existing digital 
phase-locked loops, namely to the Real-Time Combiner (RTC), the Subcarrier 
Demodulator (SD), and the Symbol Synchronization (SS) digital loops in the 
Baseband Assembly (BBA). 

II. DISCUSSION 

The phase-locked loop analyzed in this article, although totally digital, 
is in a sense a hybrid type because it operates at two extremely different 
rates. The phase detector and the integrate-and-dump circuit shown in Fig. 1 
operate at a very high rate (many millions of samples per second), while the 
loop update rate is slow (typically, one update period is longer than a 
second). Because of the very high rate at which the phase detector operates, 
the input phase, 9(t), as well as the reference phase © ( t ) are modeled as 
continuous variables. Their difference is the instantaneous loop phase error, 
namely, 


1 



INPUT PHASE 



Fig. 1. Main components of a digital phase-locked loop (DPLL) 


A 


<Kt) = e(t) - e(t) 


(i) 


In general, the phase detector is a nonlinear device: It transforms 4> C t ) to 

f (4>(t)). Two cases are most frequently considered: 


f(<j>(t)) = sin 4»(t) 


(2a) 


when tracking a sinuousoidal signal and 


f(<j>(t)) = z 


y sin(2i + 1 )<t> ( t ) 
i=0 ( 2i + l) 2 


(2b) 


when tracking a square wave. 


A typical phase error process is shown in Fig. 2. In the time interval 
(t^, t^ + ^), both <t> ( t ) and f(<j>(t)) are assumed to be monotonic functions of time. 
Because of this assumption, there is a particular instant, t^ + At^, in 
the interval (t^» t^ + ^) when the following condition is met: 


k+1 

f (4>( t k + At k ) ) = ^ / f(4>(t))dt 


(3a) 


T = t - t 
k+1 k 


2 







Figure 2. Input/output phase processes 

i.e., when f(<|>(t)) passes through its average value. In our analysis, T is 
the loop update time which is assumed to be a constant. Solving Eq. (3a), we 
define the new variable 



k+1 

£ X f(4>(t))dt 

c k 


(3b) 


We now replace the actual phase error process, <{>(t), by a staircase process 
4> k which satisfies Eq. (3b) and which is constant during one loop update 
interval. With the above substitution, Eq. (1) now becomes 


- < * >(t k ) + it k ) 


A 



( 4 ) 


3 


where 9^ is the value of 9(t) at time t^ + At^ and 9^ is the solution of 
the difference equation to be developed. 


In Fig. 1, the integrate-and-dump circuit integrates the output of the 
phase detector over an update interval producing e^, which we call the 
error signal. Because of the additive noise in the input signal, is a 
random variable. The input-output relation of the integrate-and-dump device 
can be expressed by the following equation: 

E t - V ( V + \ (5) 


where is the gain of the integrate-and-dump device defined as follows: 


G ^ 

Q " f(4>J 


( 6 ) 


The numerator, E{e^}, is the expected value of e^, and (J>^ is defined by 
Eq. (3b). The noise sample n fc is modeled as a sequence of independent 
Gaussian random variables with zero mean and variance a 2 . Table 1 gives 
expressions of G^ and o^ for various digital phase-locked loops in use at 
the Jet Propulsion Laboratory (JPL). 


In order to find the pdf of the phase error at the loop update instants, 
we need the set of difference equations that describes the operation of the 
loop. In general, between the integrate-and-dump device where the error 
signal, e, , is sampled and the digitally controlled oscillator (DCO) or 

A 

some other control device which produces the phase estimate 9(t), there is a 
digital filter which can be described by the following difference equation: 

y fc D(z) = e k N(z) (7) 

where 

D(z ) = 1 - a^z - * - ... - a n z -n , (8) 

N(z) = bp + b^z - ^ + ... + b ro z -m , m < n (9) 


4 




5 



Equation (7) 


and z is the unit delay operator defined as z 
can also be written as 


N(z ) A w v 

y k = DM C k - FU)£ 1 


( 10 ) 


where F(z) is defined as the loop filter transfer function. 


If D(z) has N integrators, then Eq. (8) can be written as 


D(z ) = (1 - z 1 ) N (1 - a. 'z 1 - ... - a 'z n+N ) (11) 

1 n-N 

The mathematical form of the difference equation for the control device in the 
feedback path will depend on the particular loop we are dealing with. If, for 
example, the loop tracks the epoch or the relative delay of a signal, as is 
the case with the RTC, then in the feedback path we have a summer , and y^ in 
Fig. 3 will be the increment in epoch (or in relative time delay). For this 
case, the input-output relationship of the control device will be described by 
the equation 


k-1 



i=-°° 



Cj^(z) 

=r y k 

(1 - z ) 


( 12 ) 


On the other hand, if the loop tracks the phase of the incoming signal with a 
DCO (carrier or subcarrier tracking loop), then y^ of Fig. 2 will be the 
frequency of the DCO. Since the phase is the integral of frequency, we can 
approximate the integral operation by the bilinear transformation (or the trap- 
ezoidal rule), namely. 


A 


e 


k 


T (1 + z 1 ) z— 1 
2 (1 - z _1 ) 


C 2 (z) 

(1 - z *) 


(13) 


Combining Eqs. (5), (10) and (12) or (13) in Eq. (4), the difference equation 
which describes the operation of the loop at the instants t^ + At is ob- 
tained, namely. 


6 




Fig. 3. Equivalent block diagram of a DPLL 


<t> k = - C^(z)F(z)(G^f (<j> k ) + n k )/(l “ z 1 )» j = 1 or 2 (14) 

The input phase, 0(t), at discrete time instants t k = kT can be modeled as 
being of the form 

0 k = e Q + 0^kT + 0 2 (kT) 2 + ... + 0 M (kT) M (15) 

where 0^, i=0, 1, 2, ...,M are constants. 

The following fact will be useful in our analysis: 

1 0 for M < N 

for M = N (16) 

unbounded otherwise 

After these preliminary steps, we are finally ready to write the overall 
nonlinear stochastic difference equation for the digital phase-locked loop. 
Using Eq. (10) in Eq. (14) and rearranging terms, we obtain 


7 






( 17 ) 


with 


<j> k (l - z 1 )D(z) = 0 k (l - z 1 )D(z) 

- C.(z)N(z)(G Q fU k ) + i^) 


j 


1 for an epoch tracking loop and 

2 for a carrier tracking loop 


Assuming that D(z) has N integrators as in Eq. (11) and that © k is of the 
form given by Eq. (15) with M = N + 1 (finite steady state error), the above 
difference equation becomes 


where 


and 


n+l n-N K 

I c(i)<j> z = u l B(i)z - l )z 


i=0 


i=0 


i=l 


K 

- I d(i)n v z 
i=l 


-l 


(18) 


K = 


m + 1 for the epoch tracking loop 
m + 2 for the carrier tracking loop 


(19) 


u £ M!* m t” (20) 

Again m and n are the number of delays in the loop filter as defined by Eqs. 

(8) and (9), and M is the highest component of (f> k tracked by the loop with 
finite steady state error. Note that c(0) = B(0) = 1. 

The above stochastic difference equation is not Markov for n larger than 
zero. However, it can easily be reduced to a set of first-order stochastic 
difference equations, resulting in a vector Markov difference equation. This 
reduction can be accomplished by several methods. Expanding the above 
difference equation and regrouping terms, we obtain 


8 


<t>k = “ 


+ z ^[J3(l)u - c(l)4> k - y( 1 )f (4> k ) ~ d(l)n k ] 

+ z _2 [J3(2)u - c( 2)<|> k - y(2)f(4> k ) - d(2)n k J + ... 

+ z n ^[13(n + l)u - c(n + 1 )4>j c - y( n + l)f($ k ) - d(n + Dn^] (21) 

(Note that I3(i) are zero for i larger than n - N, and y(i) and d(i) are zero 

for i larger than K. ) 


Representing Eq. (21) in the First Canonical Form, the flow diagram of 
Fig. 4 is obtained with the following output equation 


*k = x k (1) + u 


( 22 ) 


and the nonlinear state equation 


' x k+l(D 

X k+ i(2) 

• 

II 

-c(l) 1 0 ... 0 

-c(2) 0 1 ... 0 

0 0 ... 0 
• • • • 


x k (l) 

X k (2) 

• 

• 

X k+ i(n) 
_X k+1 (n + 1)_ 


• • • • 

-c(n) 00 1 

-c(n + 1) 0 0 0 

— J 


• 

X k (n) 

_X k (n + 1)_ 


Y ( 1 ) 

y(2) 


f(X k (l) + u) 


y(n) 

y(n + 1) 


"J3(D - c ( 1 ) 


d( 1 ) 

B(2) - c ( 2 ) 


d ( 2 ) 

• 

u - 

- 

B(n) - c(n) 


d(n) 

_ 0 - c(n + 1). 


. 0 


(23a) 


9 



n k (NOISE) 




Fig. 4. First canonical form representation of the DPLL difference equation 
In vector notation Eq. (23a) can be written as 


L = C)^. - rf(J^(l) + u) + Bu - Dn fc (23b) 

The above state equation simplifies considerably when the phase detector can 
be modeled as a linear device with f^^) = 1° this case, Eq. (18) 

reduces to 


n+1 . n-N 1 K 

I a(i)t z~ = u l fl(i)z _i - l d(i)n z (24a) 

i=0 i=0 i=l 


where 


ct( i ) = c(i) - y(i) 


(24b) 







' x k+l(D 


— ot ( 1 ) 1 0 ... 0 


x k (D 

X k (2) 

x k+l^ 2 ) 


— <x(2) 0 T IT .T 0“ 

0 0 ... 0 


X k+i (n) 


-a(n) 00 1 


X k (n) 

_X k+ i(n + 1)_ 


_-<x(n + 1)00 0 

r nr 

1 

_X k (n + 1)_ 



0(1) - a(l ) 
0(2) - a(2 ) 


d(l) 

d(2) 

+ 

0(n) - a(n) 

u - 

d(n) 


0 - a(n + 1 ) . 


. o . 


n k 


(25a) 


Using vector notation we can write Eq. (25a) 


X. , = AX, + Bu - Dn, 
— k+1 -k - - k 


(25b) 


where A is an (n + 1) x (n + 1) matrix. The values of the constants {a^}, 
{13^}, and {d^} corresponding to the RTC and the subcarrier loop are given 
in Table 2. 


A. Chapman-Kolmogorov Equation 


The vector Markov process for a nonlinear phase detector is given by 
Eq. (23) and for a linear phase detector by Eq. (26). The conditional proba- 
bility density function (pdf) of X^ + ^, conditioned on the initial condition 
Xq, satisfies the Chapman-Kolmogorov (C-K) equation 


- 


00 CD 

/ .... / % +l = SlSk 

_CO —CO 

x P k^ = -l-0* dw l *** dW L 


w) 


where 


Xq = initial value of X 
Pk^ = x|x Q ) = Pdf of X^ given x Q 
^k^^k+1 = = = transition pdf of X^ + ^ given X^ = w 


(26) 


LI 


Table 2. Value of constants for the linear state Eq. (24) 


Constants 

RTC 

Subcarrier and Symbol Synchronization 
Loops 

«1 

AGq - 1 - B 

TAGq/2 - C - 2 

a 2 

B 

2C + 1 + TGq(A + B)/2 

a 3 

— 

2D - C + TGqB/2 

a 4 

— 

-D 

(3 0 

1 

1 

h 

-B 

-C 

A 2 

0 

-D 

b 3 

— 

0 


— 

0 

<*1 

-A 

-TA/2 

d 2 

0 

-T(A + B)/2 

d 3 

— 

-TB/2 

d 4 

— 

0 
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Now. the- transition_pdf_.q^ (_? !_• )_can b_e_ written _as _ 

= q k lx k+i (1) = x(1) * \+i (2) = x(2) * •••• \+i a) = x(L) |4 = ^ 

- q k [X k+i (L) = X (L)|x k+1 (L - 1) = x(L - 1) X k+ i(i) = x(l) , X k = w] 

x q k [x k+i (L " 15 = x(L " 1 ) l x k+i (L - 2 ) “ *a - 2) x k+ i (1) = x(1) > x k - 


x q k [x k+i (1) = x(1) Uk = 
where L = n + 1 . 

In words, Eq. (27) expresses the total transition pdf as a product of 
transition pdf's at the input to each delay of Fig. 4. If we refer to that 
figure and assume Gaussian distribution of the noise samples n fc , these 
transition pdf's will have Gaussian distributions for d(i) ^ 0 and will be 
delta functions for d(i) =0, i = 1, n + 1. The pdf's with Gaussian 

distribution will have the following conditional mean and variance 

E ^ X k+l (i) l X k+l (i ” L) = x(i ' 1)s •**» x k+l (1) = x(1) ’ \ = = M i 

= -c(i)w(l) -y(i)f(w(l) + u) + w(i + 1) +(Q(i )-c(i ) )u 

(28) 

2 2 2 

of = d(i) a (29) 

l n 

and the delta functions will be of the form 


5I W i) - ''i 1 


(30) 



Note in Eqs . (28), (29), and (30) that the transition pdf (•(•) 

(Eq. [27]) will be independent of k. This should be expected since we are 
analyzing a time-invariant system. For simplicity of notation, we shall define 

q k (j) = q k [X k+l (j)|X k + l (j-1) = *••• \ + l (1) = x(1) ’ \ ^ 

j = L, L-l , ..., 1 (31) 

The transition pdf of Eq. (27) now can be written as 

\ (i> (32) 
3=1 

The C-K equation for an arbitrary nonlinearity f(») is given by Eq. (26). 

When f(«) is periodic, as in Eq. (2), then the pdf of the unrestricted phase 
error process p^ + ^(X^ + ^ = x|X = w) will go to zero as k ■* °°, which is not 
particularly useful for loop analysis. However, the modulo 2ir reduction of 
the random variable is of interest in studying the behavior of the loop in 

the locked condition. We define a new random variable ’ ( j ) e(-ir, ir) 
which is obtained from the old random variable X^j) by taking 

^'(j) = ^(j) mod 2ir, j = 1,2 L (33) 

Using Eqs. (26) and (32), the corresponding C-K equation that applies to the 
mod 2ir process {X’} is modified to be 

p k+l ( W = S **• S " q k (j)p k ( V = w|x 0 )dw l--- dw L (34) 

-ir -ir 3=1 

where all the variables are now e(-ir, ir). The transition pdf's, 9 k (j) 
of Eq. (31), will have to be modified now as follows: 


14 





m=-°° 

= x(j-l), ..., X^' = x(l), X^ ' = w] (35) 

In solving the C-K equation as expressed by Eq. (34) in the steady state, 
i.e., k ■* 00 , we know that p^ + ^( # ) on the left hand side will be equal to 
p^(») on the right hand side. The dependence of distribution p^(*) on 
Xq, the initial value of the state vector, can be dropped for large k. 

In order to conclude our discussion, we will show a method of solving the 
C-K equation for the first-order loop, i.e., when L = (n + 1) = 1 (no loop 
filter). In this case the transition pdf will have only one term with 
Gaussian distribution, namely. 


V *"> 


00 



exp 


( (X* +2irm) - 


( — cw - yf (w+u)-cu)}‘ 
2 o 2 


(36) 


Note that we have dropped the index 1 from all the variables and constants in 
order to simplify our notation. 

The C-K equation for this first-order loop in the steady state will be 
-tr 00 

p(X' = x) = , = S I exp 

J 2-rro -ir m=-°° 

X' e (-T , n ) (37) 

We now approximate integration with summation. The continuous variables X' and 

w are converted to discrete variables taking values x. and w. in the 

ii. 

interval (-tr, ir), namely. 


{(X'+2mn) - ( — cw - yf (w+u)-cu)}‘ 
2 a 2 


p(X = w) dw, 


x^ = i(tr/N), w^ = i(ir/N), i = -N,-N+l, ..., N 


(39) 


where N is an integer large enough to make the above approximation valid. 
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With these changes, the discrete C-K equation becomes 


p(X' 


N 

= x i) = l q^pOt' 

j=-N J 


w. ) 
J 


where 


( 40 ) 


ij 


A 


exp 


2iro m=- 




{(x.+2irm) - (-cw. - yf (w ,+u)-cu)}‘ 
2o 2 




(41) 


Let Q be the (2N + 1) x (2N + 1) matrix with elements q„ defined by Eq. 
(41), and let P be the (2N +1) vector with elements x^ and w^. Then the 
discrete C-K equation in vector notation will be 


P = QP 


(42) 


Using the definition of eigenvalues, we recognize that Eq. (42) can be 
written as 


(XI - Q)P = 0 (43) 

where X is an eigenvalue of Q and I is a (2N+l)x(2N+ 1) identity matrix. 
The C-K discrete equation will have a solution if one of the eigenvalues of Q 
will be 1. In this case, the desired discrete pdf of the phase error P will 
be the eigenvector corresponding to X = 1. 

B. A Numerical Example 

The modulo 2ir steady-state pdf of the phase error process was 
numerically computed for a first order loop having a summer as the control 
device and a filter with F(z) = A. It was assumed that the loop tracks a 
square wave which results in the sawtooth nonlinearity given by Eq. (2). 

These assumptions result in the following nonlinear stochastic difference 
equation: 


♦k - s k - z ~ 1a <V <v - V /(1 - 


(44) 
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-The elements q.j-of the -matrix Q -defined- by Eq. -(42) are -normalized- as 
follows 


q ij 


q„ir/(2N + 1) for i = 1 and (2N + 1) 
q„2ir/(2N + 1) otherwise 


(45) 


where q„ are computed using Eq. (41). This normalization makes the sum of 
the elements q„' along each column of Q' equal to 1.0 which in turn makes 
the largest eigenvalue of Q', , equal to 1.0 and the desired, discrete 

pdf, P, is obtained. In our numerical evaluation of P, the following assump- 
tions were made: 


0.01 < A < 1.0 

u = 0.0 (no loop stress) 

0.01 < a 2 < 1.0 

Figures 5, 6, and 7 summarize the numerical results. In these figures the 
pdfs of the modulo 2ir process of the phase error are compared with a 
Gaussian pdf of the same mean and variance values. 

III. Conclusions 

Using a method based on transforming high-order nonlinear stochastic 
difference equations to a set of first-order nonlinear difference equations, 
the steady-state probability density function of a digital phase-locked loop 
with integrate-and-dump circuit following the phase detector has been derived. 
This work ties together in a unified fashion, the previous analyses of specific 
DSN Baseband Assembly (BBA) loops such as the Real-Time Combiner (RTC), Sub- 
carrier Demodulation Assembly (SDA), and Symbol Synchronization Assembly (SSA). 
A numerical example for a first-order loop having a summer as a control device 
and a constant loop filter was given to illustrate the simplicity of the appli- 
cation of the theory. It is believed that the results derived here are suffi- 
ciently general so as to apply to a host of other loops of this generic type. 
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Fig. 6. Modulo 2 t r steady-state pdf of the phase error. First order loop. 
Loop gain GqA = 0.1. 
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